df_uk_covid$time %>% summary()
Min. 1st Qu. Median Mean 3rd Qu. Max.
1.00 8.75 16.50 16.50 24.25 32.00
df_uk_ctrl_nuts <- read_csv("controls_UK_nuts3.csv")
Parsed with column specification:
cols(
nuts3 = [31mcol_character()[39m,
nuts3_name = [31mcol_character()[39m,
airport_dist = [32mcol_double()[39m,
males = [32mcol_double()[39m,
popdens = [32mcol_double()[39m,
manufacturing = [32mcol_double()[39m,
tourism = [32mcol_double()[39m,
health = [32mcol_double()[39m,
academic = [32mcol_double()[39m,
medinc = [32mcol_double()[39m,
medage = [32mcol_double()[39m,
conservative = [32mcol_double()[39m
)
df_uk_ctrl_nuts <- df_uk_ctrl_nuts %>% select(-nuts3_name)
df_uk_ctrl_nuts
df_uk_ctrl_ut <- read_csv("controls_UK_ut.csv")
Parsed with column specification:
cols(
ut_area = [31mcol_character()[39m,
ut_name = [31mcol_character()[39m,
airport_dist = [32mcol_double()[39m,
males = [32mcol_double()[39m,
popdens = [32mcol_double()[39m,
manufacturing = [32mcol_double()[39m,
tourism = [32mcol_double()[39m,
health = [32mcol_double()[39m,
academic = [32mcol_double()[39m,
medinc = [32mcol_double()[39m,
medage = [32mcol_double()[39m,
conservative = [32mcol_double()[39m
)
df_uk_ctrl_ut <- df_uk_ctrl_ut %>% select(-ut_name)
df_uk_ctrl_ut
NA
NA
df_uk <- df_uk %>% plyr::join(df_uk_ctrl_ut, by='ut_area')
df_uk_socdist <- df_uk_socdist %>% plyr::join(df_uk_ctrl_nuts, by='nuts3')
nuts_london_inner <- c('UKI31','UKI32','UKI33','UKI34','UKI41',
'UKI42','UKI43','UKI44','UKI45')
nuts_london_outer <- c('UKI51','UKI52','UKI53','UKI54','UKI61',
'UKI62','UKI63','UKI71','UKI72','UKI73',
'UKI74','UKI75')
ut_london_inner <- c('E09000007','E09000001','E09000033','E09000013',
'E09000020','E09000032','E09000025','E09000012',
'E09000030','E09000014','E09000019','E09000023',
'E09000028','E09000022')
ut_london_outer <- c('E09000011','E09000004','E09000016','E09000002',
'E09000031','E09000026','E09000010','E09000006',
'E09000008','E09000029','E09000021','E09000024',
'E09000003','E09000005','E09000009','E09000017',
'E09000015','E09000018','E09000027')
df_uk = df_uk %>%
mutate(london = ifelse(ut_area %in% ut_london_inner, 'london_inner',
ifelse(ut_area %in% ut_london_outer, 'london_outer',
'country'))) %>%
mutate(london = as.factor(london))
df_uk_socdist = df_uk_socdist %>%
mutate(london = ifelse(nuts3 %in% nuts_london_inner, 'london_inner',
ifelse(nuts3 %in% nuts_london_outer, 'london_outer',
'country'))) %>%
mutate(london = as.factor(london))
df_uk %>% ggplot(aes(x=time, y=rate_day)) +
geom_point(aes(col=ut_area, size=popdens)) +
geom_smooth(method="loess", se=T) +
theme(legend.position="none") +
ggtitle("Overall prevalence over time")
pers <- c('pers_o', 'pers_c', 'pers_e', 'pers_a', 'pers_n')
for (i in pers){
gg <- df_uk %>% mutate(prev_tail = cut(.[[i]],
breaks = c(-Inf, quantile(.[[i]], 0.2), quantile(.[[i]], 0.8), Inf),
labels = c('lower tail', 'center', 'upper tail'))) %>%
filter(prev_tail != 'center') %>%
ggplot(aes(x=time, y=rate_day)) +
geom_point(aes(col=ut_area, size=popdens)) +
geom_smooth(method="loess", se=T) +
facet_wrap(~prev_tail) +
theme(legend.position="none") +
ggtitle(i)
print(gg)
}
df_uk %>% ggplot(aes(x=time, y=rate_day)) +
geom_point(aes(col=ut_area, size=popdens)) +
geom_smooth(method="loess", se=T) +
theme(legend.position="none") +
facet_wrap(~london) +
ggtitle("Overall prevalence over time")
NA
NA
NA
df_uk_socdist %>% ggplot(aes(x=time, y=socdist_single_tile)) +
geom_point(aes(col=nuts3, size=popdens)) +
geom_smooth(method="loess", se=T) +
theme(legend.position="none") +
facet_wrap(~london) +
ggtitle("Overall social distancing over time")
weekend <- c(6, 7, 13, 14, 20, 21)
df_uk_loess <- df_uk_socdist %>% filter(!time %in% weekend) %>%
split(.$nuts3) %>%
map(~ loess(socdist_single_tile ~ time, data = .)) %>%
map(predict, 1:23) %>%
bind_rows() %>%
gather(key = 'nuts3', value = 'loess') %>%
group_by(nuts3) %>%
mutate(time = row_number())
df_uk_socdist <- df_uk_socdist %>% merge(df_uk_loess, by=c('nuts3', 'time')) %>%
mutate(socdist_single_tile_clean = ifelse(time %in% weekend, loess,
socdist_single_tile)) %>%
arrange(nuts3, time)
df_uk_socdist %>% ggplot(aes(x=time, y=loess, group=nuts3)) +
geom_line()
df_uk_socdist
NA
NA
df_uk_socdist %>% ggplot(aes(x=time, y=socdist_single_tile_clean)) +
geom_point(aes(col=nuts3, size=popdens)) +
geom_smooth(method="loess", se=T) +
theme(legend.position="none") +
ggtitle("Overall social distancing over time")
pers <- c('pers_o', 'pers_c', 'pers_e', 'pers_a', 'pers_n')
for (i in pers){
gg <- df_uk_socdist %>% mutate(socdist_tail = cut(.[[i]],
breaks = c(-Inf, quantile(.[[i]], 0.2), quantile(.[[i]], 0.8), Inf),
labels = c('lower tail', 'center', 'upper tail'))) %>%
filter(socdist_tail != 'center') %>%
ggplot(aes(x=time, y=socdist_single_tile_clean)) +
geom_point(aes(col=nuts3, size=popdens)) +
geom_smooth(method="loess", se=T) +
facet_wrap(~socdist_tail) +
theme(legend.position="none") +
ggtitle(i)
print(gg)
}
df_uk_socdist <- df_uk_socdist %>% mutate(socdist_single_tile = socdist_single_tile_clean) %>%
select(-loess, -socdist_single_tile_clean)
df_uk %>% group_by(ut_area) %>%
summarize_if(is.numeric, mean, na.rm=T) %>%
select(-ut_area, -time) %>%
cor(use = 'pairwise.complete') %>% round(3)
pers_o pers_e pers_a pers_n pers_c frequ rate_day airport_dist males popdens
pers_o 1.000 0.734 -0.604 -0.143 -0.636 0.011 0.647 -0.247 0.437 0.806
pers_e 0.734 1.000 -0.412 -0.497 -0.293 0.071 0.559 -0.262 0.197 0.617
pers_a -0.604 -0.412 1.000 -0.197 0.643 0.130 -0.584 0.285 -0.460 -0.745
pers_n -0.143 -0.497 -0.197 1.000 -0.353 -0.243 -0.108 0.059 0.096 0.043
pers_c -0.636 -0.293 0.643 -0.353 1.000 0.220 -0.487 0.327 -0.524 -0.727
frequ 0.011 0.071 0.130 -0.243 0.220 1.000 -0.077 0.108 -0.124 -0.180
rate_day 0.647 0.559 -0.584 -0.108 -0.487 -0.077 1.000 -0.388 0.302 0.758
airport_dist -0.247 -0.262 0.285 0.059 0.327 0.108 -0.388 1.000 -0.124 -0.381
males 0.437 0.197 -0.460 0.096 -0.524 -0.124 0.302 -0.124 1.000 0.510
popdens 0.806 0.617 -0.745 0.043 -0.727 -0.180 0.758 -0.381 0.510 1.000
manufacturing -0.532 -0.595 0.496 0.359 0.363 -0.111 -0.511 0.300 -0.259 -0.545
tourism 0.116 0.098 0.069 -0.244 0.227 0.115 -0.023 0.448 -0.094 -0.099
health 0.089 0.026 -0.001 0.138 -0.286 -0.132 -0.010 0.150 -0.012 0.112
academic 0.707 0.730 -0.502 -0.447 -0.298 0.152 0.647 -0.369 0.267 0.602
medinc 0.533 0.569 -0.533 -0.335 -0.300 0.070 0.573 -0.379 0.367 0.631
medage -0.508 -0.337 0.582 -0.167 0.780 0.144 -0.513 0.447 -0.620 -0.705
conservative -0.843 -0.767 0.511 0.366 0.485 -0.125 -0.626 0.396 -0.293 -0.686
manufacturing tourism health academic medinc medage conservative
pers_o -0.532 0.116 0.089 0.707 0.533 -0.508 -0.843
pers_e -0.595 0.098 0.026 0.730 0.569 -0.337 -0.767
pers_a 0.496 0.069 -0.001 -0.502 -0.533 0.582 0.511
pers_n 0.359 -0.244 0.138 -0.447 -0.335 -0.167 0.366
pers_c 0.363 0.227 -0.286 -0.298 -0.300 0.780 0.485
frequ -0.111 0.115 -0.132 0.152 0.070 0.144 -0.125
rate_day -0.511 -0.023 -0.010 0.647 0.573 -0.513 -0.626
airport_dist 0.300 0.448 0.150 -0.369 -0.379 0.447 0.396
males -0.259 -0.094 -0.012 0.267 0.367 -0.620 -0.293
popdens -0.545 -0.099 0.112 0.602 0.631 -0.705 -0.686
manufacturing 1.000 -0.088 -0.023 -0.714 -0.532 0.472 0.657
tourism -0.088 1.000 0.049 0.056 -0.062 0.435 -0.023
health -0.023 0.049 1.000 -0.228 -0.243 -0.153 0.043
academic -0.714 0.056 -0.228 1.000 0.733 -0.382 -0.888
medinc -0.532 -0.062 -0.243 0.733 1.000 -0.481 -0.661
medage 0.472 0.435 -0.153 -0.382 -0.481 1.000 0.486
conservative 0.657 -0.023 0.043 -0.888 -0.661 0.486 1.000
df_uk_socdist %>% group_by(nuts3) %>%
summarize_if(is.numeric, mean, na.rm=T) %>%
select(-nuts3, -time) %>%
cor(use = 'pairwise.complete') %>% round(3)
rate_day socdist_tiles socdist_single_tile pers_o pers_e pers_a pers_n pers_c frequ
rate_day 1.000 -0.695 0.487 0.635 0.562 -0.567 -0.113 -0.469 0.078
socdist_tiles -0.695 1.000 -0.619 -0.626 -0.733 0.545 0.325 0.395 -0.311
socdist_single_tile 0.487 -0.619 1.000 0.317 0.356 -0.382 0.040 -0.298 0.032
pers_o 0.635 -0.626 0.317 1.000 0.713 -0.603 -0.122 -0.651 0.183
pers_e 0.562 -0.733 0.356 0.713 1.000 -0.450 -0.456 -0.363 0.249
pers_a -0.567 0.545 -0.382 -0.603 -0.450 1.000 -0.215 0.661 -0.030
pers_n -0.113 0.325 0.040 -0.122 -0.456 -0.215 1.000 -0.338 -0.332
pers_c -0.469 0.395 -0.298 -0.651 -0.363 0.661 -0.338 1.000 0.096
frequ 0.078 -0.311 0.032 0.183 0.249 -0.030 -0.332 0.096 1.000
airport_dist -0.353 0.550 -0.357 -0.231 -0.319 0.328 0.072 0.326 -0.147
males 0.303 -0.287 0.077 0.519 0.243 -0.567 0.195 -0.632 -0.032
popdens 0.730 -0.705 0.500 0.803 0.638 -0.748 0.030 -0.737 -0.040
manufacturing -0.412 0.668 -0.438 -0.496 -0.605 0.427 0.343 0.309 -0.278
tourism -0.048 0.239 -0.110 0.079 0.005 0.142 -0.168 0.207 -0.007
health 0.002 0.140 -0.060 0.086 0.047 0.005 0.151 -0.274 -0.258
academic 0.640 -0.808 0.348 0.713 0.742 -0.491 -0.432 -0.324 0.427
medinc 0.600 -0.772 0.457 0.552 0.582 -0.620 -0.247 -0.367 0.248
medage -0.472 0.563 -0.328 -0.512 -0.402 0.631 -0.181 0.790 -0.047
conservative -0.629 0.757 -0.327 -0.837 -0.761 0.514 0.338 0.522 -0.374
airport_dist males popdens manufacturing tourism health academic medinc medage
rate_day -0.353 0.303 0.730 -0.412 -0.048 0.002 0.640 0.600 -0.472
socdist_tiles 0.550 -0.287 -0.705 0.668 0.239 0.140 -0.808 -0.772 0.563
socdist_single_tile -0.357 0.077 0.500 -0.438 -0.110 -0.060 0.348 0.457 -0.328
pers_o -0.231 0.519 0.803 -0.496 0.079 0.086 0.713 0.552 -0.512
pers_e -0.319 0.243 0.638 -0.605 0.005 0.047 0.742 0.582 -0.402
pers_a 0.328 -0.567 -0.748 0.427 0.142 0.005 -0.491 -0.620 0.631
pers_n 0.072 0.195 0.030 0.343 -0.168 0.151 -0.432 -0.247 -0.181
pers_c 0.326 -0.632 -0.737 0.309 0.207 -0.274 -0.324 -0.367 0.790
frequ -0.147 -0.032 -0.040 -0.278 -0.007 -0.258 0.427 0.248 -0.047
airport_dist 1.000 -0.155 -0.377 0.272 0.492 0.180 -0.379 -0.389 0.469
males -0.155 1.000 0.615 -0.180 -0.171 0.028 0.263 0.440 -0.672
popdens -0.377 0.615 1.000 -0.492 -0.135 0.168 0.574 0.652 -0.714
manufacturing 0.272 -0.180 -0.492 1.000 -0.051 -0.069 -0.640 -0.462 0.413
tourism 0.492 -0.171 -0.135 -0.051 1.000 0.027 0.005 -0.138 0.485
health 0.180 0.028 0.168 -0.069 0.027 1.000 -0.202 -0.261 -0.154
academic -0.379 0.263 0.574 -0.640 0.005 -0.202 1.000 0.724 -0.368
medinc -0.389 0.440 0.652 -0.462 -0.138 -0.261 0.724 1.000 -0.479
medage 0.469 -0.672 -0.714 0.413 0.485 -0.154 -0.368 -0.479 1.000
conservative 0.389 -0.350 -0.678 0.627 0.012 0.010 -0.891 -0.653 0.499
conservative
rate_day -0.629
socdist_tiles 0.757
socdist_single_tile -0.327
pers_o -0.837
pers_e -0.761
pers_a 0.514
pers_n 0.338
pers_c 0.522
frequ -0.374
airport_dist 0.389
males -0.350
popdens -0.678
manufacturing 0.627
tourism 0.012
health 0.010
academic -0.891
medinc -0.653
medage 0.499
conservative 1.000
# function calculates all relevant models
run_models <- function(y, lvl1_x, lvl2_x, lvl2_id, data, ctrls=F){
# subset data
data = data %>%
dplyr::select(all_of(y), all_of(lvl1_x), all_of(lvl2_x), all_of(lvl2_id),
popdens, rate_day, all_of(y))
data = data %>%
dplyr::rename(y = all_of(y),
lvl1_x = all_of(lvl1_x),
lvl2_x = all_of(lvl2_x),
lvl2_id = all_of(lvl2_id)
)
# configure optimization procedure
ctrl_config <- lmeControl(opt = 'optim', maxIter = 100, msMaxIter = 100)
# baseline
baseline <- lme(fixed = y ~ 1, random = ~ 1 | lvl2_id,
data = data,
correlation = corAR1(),
control = ctrl_config,
method = 'ML')
# random intercept fixed slope
random_intercept <- lme(fixed = y ~ lvl1_x + lvl2_x,
random = ~ 1 | lvl2_id,
data = data,
correlation = corAR1(),
control = ctrl_config,
method = 'ML')
# random intercept random slope
random_slope <- lme(fixed = y ~ lvl1_x + lvl2_x,
random = ~ lvl1_x | lvl2_id,
data = data,
correlation = corAR1(),
control = ctrl_config,
method = 'ML')
# cross level interaction
interaction <- lme(fixed = y ~ lvl1_x * lvl2_x,
random = ~ lvl1_x | lvl2_id,
data = data,
correlation = corAR1(),
control = ctrl_config,
method = 'ML')
# create list with results
results <- list('baseline' = baseline,
"random_intercept" = random_intercept,
"random_slope" = random_slope,
"interaction" = interaction)
if (ctrls == 'dem' | ctrls == 'prev'){
# random intercept random slope
random_slope_ctrl_dem <- lme(fixed = y ~ lvl1_x + lvl2_x + popdens,
random = ~ lvl1_x | lvl2_id,
data = data,
correlation = corAR1(),
control = ctrl_config,
method = 'ML')
# cross level interaction
interaction_ctrl_main_dem <- lme(fixed = y ~ lvl1_x * lvl2_x + popdens,
random = ~ lvl1_x | lvl2_id,
data = data,
correlation = corAR1(),
control = ctrl_config,
method = 'ML')
# cross level interaction
interaction_ctrl_int_dem <- lme(fixed = y ~ lvl1_x * lvl2_x + lvl1_x * popdens,
random = ~ lvl1_x | lvl2_id,
data = data,
correlation = corAR1(),
control = ctrl_config,
method = 'ML')
# create list with results
results <- list('baseline' = baseline,
"random_intercept" = random_intercept,
"random_slope" = random_slope,
"interaction" = interaction,
"random_slope_ctrl_dem" = random_slope_ctrl_dem,
"interaction_ctrl_main_dem" = interaction_ctrl_main_dem,
"interaction_ctrl_int_dem" = interaction_ctrl_int_dem)
}
if (ctrls == 'prev'){
# random intercept random slope
random_slope_ctrl_prev <- lme(fixed = y ~ lvl1_x + lvl2_x + popdens + rate_day,
random = ~ lvl1_x + rate_day | lvl2_id,
data = data,
correlation = corAR1(),
control = ctrl_config,
method = 'ML')
# cross level interaction
interaction_ctrl_main_prev <- lme(fixed = y ~ lvl1_x * lvl2_x + popdens + rate_day,
random = ~ lvl1_x | lvl2_id,
data = data,
correlation = corAR1(),
control = ctrl_config,
method = 'ML')
# cross level interaction
interaction_ctrl_int_prev<- lme(fixed = y ~ lvl1_x * lvl2_x + lvl1_x * popdens + rate_day,
random = ~ lvl1_x + rate_day | lvl2_id,
data = data,
correlation = corAR1(),
control = ctrl_config,
method = 'ML')
# create list with results
results <- list('baseline' = baseline,
"random_intercept" = random_intercept,
"random_slope" = random_slope,
"interaction" = interaction,
"random_slope_ctrl_dem" = random_slope_ctrl_dem,
"interaction_ctrl_main_dem" = interaction_ctrl_main_dem,
"interaction_ctrl_int_dem" = interaction_ctrl_int_dem,
"random_slope_ctrl_prev" = random_slope_ctrl_prev,
"interaction_ctrl_main_prev" = interaction_ctrl_main_prev,
"interaction_ctrl_int_prev" = interaction_ctrl_int_prev)
}
if(ctrls == 'exp'){
# random intercept random slope
random_slope_exp <- lme(fixed = y ~ (lvl1_x + I(lvl1_x^2)) + lvl2_x,
random = ~ (lvl1_x + I(lvl1_x^2)) | lvl2_id,
data = data,
correlation = corAR1(),
control = ctrl_config,
method = 'ML')
# cross level interaction
interaction_exp <- lme(fixed = y ~ (lvl1_x + I(lvl1_x^2)) * lvl2_x,
random = ~ (lvl1_x + I(lvl1_x^2)) | lvl2_id,
data = data,
correlation = corAR1(),
control = ctrl_config,
method = 'ML')
# create list with results
results <- list('baseline' = baseline,
"random_intercept" = random_intercept,
"random_slope" = random_slope,
"interaction" = interaction,
"random_slope_exp" = random_slope_exp,
"interaction_exp" = interaction_exp)
}
return(results)
}
# extracts table with coefficients and tests statistics
extract_results <- function(models) {
models_summary <- models %>%
map(summary) %>%
map("tTable") %>%
map(as.data.frame) %>%
map(round, 10)
# %>% map(~ .[str_detect(rownames(.), 'Inter|lvl'),])
return(models_summary)
}
# calculates comparison of all models in model list
compare_models <- function(models) {
mdl_names <- models %>% names()
str = ''
for (i in mdl_names){
mdl_str <- paste('models$', i, sep = '')
if(i == 'baseline'){
str <- mdl_str
}else{
str <- paste(str, mdl_str, sep=', ')
}
}
anova_str <- paste0('anova(', str, ')')
mdl_comp <- eval(parse(text=anova_str))
rownames(mdl_comp) = mdl_names
return(mdl_comp)
}
# df_uk <- df_uk %>% filter(london == 'country')
# df_uk_socdist <- df_uk_socdist %>% filter(london == 'country')
lvl2_scaled_ut <- df_uk %>%
dplyr::select(-time, -frequ, -rate_day, -london) %>%
distinct() %>%
mutate_at(vars(-ut_area), scale)
lvl1_scaled_ut <- df_uk %>% select(ut_area, time, rate_day) %>%
mutate_at(vars(-ut_area, -time), scale)
df_uk_scaled <- plyr::join(lvl1_scaled_ut, lvl2_scaled_ut, by = 'ut_area')
head(df_uk_scaled)
lvl2_scaled_nuts <- df_uk_socdist %>%
dplyr::select(-time, -date, -frequ, -london,
-socdist_tiles, -socdist_single_tile, -rate_day) %>%
distinct() %>%
mutate_at(vars(-nuts3), scale)
lvl1_scaled_nuts <- df_uk_socdist %>% select(nuts3, time, socdist_single_tile, rate_day) %>%
mutate_at(vars(-nuts3, -time), scale)
df_uk_socdist_scaled <- plyr::join(lvl1_scaled_nuts, lvl2_scaled_nuts, by = 'nuts3')
head(df_uk_socdist_scaled)
NA
models_o_covid <-run_models(y = 'rate_day',
lvl1_x = 'time',
lvl2_x = 'pers_o',
lvl2_id = 'ut_area',
data = df_uk_scaled,
ctrls = 'dem')
extract_results(models_o_covid)
$baseline
$random_intercept
$random_slope
$interaction
$random_slope_ctrl_dem
$interaction_ctrl_main_dem
$interaction_ctrl_int_dem
compare_models(models_o_covid)
NA
models_c_covid <-run_models(y = 'rate_day',
lvl1_x = 'time',
lvl2_x = 'pers_c',
lvl2_id = 'ut_area',
data = df_uk_scaled,
ctrls = 'dem')
extract_results(models_c_covid)
$baseline
$random_intercept
$random_slope
$interaction
$random_slope_ctrl_dem
$interaction_ctrl_main_dem
$interaction_ctrl_int_dem
compare_models(models_c_covid)
NA
NA
models_e_covid <-run_models(y = 'rate_day',
lvl1_x = 'time',
lvl2_x = 'pers_e',
lvl2_id = 'ut_area',
data = df_uk_scaled,
ctrls = 'dem')
extract_results(models_e_covid)
$baseline
$random_intercept
$random_slope
$interaction
$random_slope_ctrl_dem
$interaction_ctrl_main_dem
$interaction_ctrl_int_dem
compare_models(models_e_covid)
NA
NA
models_a_covid <-run_models(y = 'rate_day',
lvl1_x = 'time',
lvl2_x = 'pers_a',
lvl2_id = 'ut_area',
data = df_uk_scaled,
ctrls = 'dem')
extract_results(models_a_covid)
$baseline
$random_intercept
$random_slope
$interaction
$random_slope_ctrl_dem
$interaction_ctrl_main_dem
$interaction_ctrl_int_dem
compare_models(models_a_covid)
NA
NA
models_n_covid <-run_models(y = 'rate_day',
lvl1_x = 'time',
lvl2_x = 'pers_n',
lvl2_id = 'ut_area',
data = df_uk_scaled,
ctrls = 'dem')
extract_results(models_n_covid)
$baseline
$random_intercept
$random_slope
$interaction
$random_slope_ctrl_dem
$interaction_ctrl_main_dem
$interaction_ctrl_int_dem
compare_models(models_n_covid)
NA
NA
models_o_covid_exp <-run_models(y = 'rate_day',
lvl1_x = 'time',
lvl2_x = 'pers_o',
lvl2_id = 'ut_area',
data = df_uk_scaled,
ctrls = 'exp')
extract_results(models_o_covid_exp)
$baseline
$random_intercept
$random_slope
$interaction
$random_slope_exp
$interaction_exp
compare_models(models_o_covid_exp)
NA
models_c_covid_exp <-run_models(y = 'rate_day',
lvl1_x = 'time',
lvl2_x = 'pers_c',
lvl2_id = 'ut_area',
data = df_uk_scaled,
ctrls = 'exp')
extract_results(models_c_covid_exp)
$baseline
$random_intercept
$random_slope
$interaction
$random_slope_exp
$interaction_exp
compare_models(models_c_covid_exp)
NA
models_e_covid_exp <-run_models(y = 'rate_day',
lvl1_x = 'time',
lvl2_x = 'pers_e',
lvl2_id = 'ut_area',
data = df_uk_scaled,
ctrls = 'exp')
extract_results(models_e_covid_exp)
$baseline
$random_intercept
$random_slope
$interaction
$random_slope_exp
$interaction_exp
compare_models(models_e_covid_exp)
NA
models_a_covid_exp <-run_models(y = 'rate_day',
lvl1_x = 'time',
lvl2_x = 'pers_a',
lvl2_id = 'ut_area',
data = df_uk_scaled,
ctrls = 'exp')
extract_results(models_a_covid_exp)
$baseline
$random_intercept
$random_slope
$interaction
$random_slope_exp
$interaction_exp
compare_models(models_a_covid_exp)
NA
models_n_covid_exp <-run_models(y = 'rate_day',
lvl1_x = 'time',
lvl2_x = 'pers_n',
lvl2_id = 'ut_area',
data = df_uk_scaled,
ctrls = 'exp')
extract_results(models_n_covid_exp)
$baseline
$random_intercept
$random_slope
$interaction
$random_slope_exp
$interaction_exp
compare_models(models_n_covid_exp)
NA
summary_table <- function(models, dv_name, prev=F){
temp_df_ctrl_main <- NULL
temp_df_ctrl_int <- NULL
temp_df_ctrl_int_prev <- NULL
for (i in models){
results <- i %>% extract_results()
results_ctrl_main <- results$interaction_ctrl_main_dem['lvl1_x:lvl2_x',]
temp_df_ctrl_main <- temp_df_ctrl_main %>% rbind(results_ctrl_main)
results_ctrl_int <- results$interaction_ctrl_int_dem['lvl1_x:lvl2_x',]
temp_df_ctrl_int <- temp_df_ctrl_int %>% rbind(results_ctrl_int)
if(prev){
results_ctrl_int_prev <- results$interaction_ctrl_int_prev['lvl1_x:lvl2_x',]
temp_df_ctrl_int_prev <- temp_df_ctrl_int_prev %>% rbind(results_ctrl_int_prev)
}
}
names_ctrl_main <- paste0(dv_name, '~', c('o', 'c', 'e', 'a', 'n'), '*time', '_crtl_popdens')
rownames(temp_df_ctrl_main) <- names_ctrl_main
names_ctrl_int <- paste0(dv_name, '~', c('o', 'c', 'e', 'a', 'n'), '*time', '_crtl_popdens*time')
rownames(temp_df_ctrl_int) <- names_ctrl_int
if(prev){
names_ctrl_int_prev <- paste0(dv_name, '~', c('o', 'c', 'e', 'a', 'n'), '*time', '_crtl_popdens*time_prev')
rownames(temp_df_ctrl_int_prev) <- names_ctrl_int_prev
sum_tab <- rbind(temp_df_ctrl_main, temp_df_ctrl_int, temp_df_ctrl_int_prev) %>% round(4)
}else{
sum_tab <- rbind(temp_df_ctrl_main, temp_df_ctrl_int) %>% round(4)
}
return(sum_tab)
}
# prevalence
models_prev <- list(models_o_covid,
models_c_covid,
models_e_covid,
models_a_covid,
models_n_covid)
sum_tab_prev <- summary_table(models_prev, dv_name = 'prev')
write.table(sum_tab_prev, quote=F)
Value Std.Error DF t-value p-value
prev~o*time_crtl_popdens 0.0448 0.0025 3127 17.8777 0
prev~c*time_crtl_popdens -0.0386 0.0025 3127 -15.1522 0
prev~e*time_crtl_popdens 0.0372 0.0026 3127 14.5749 0
prev~a*time_crtl_popdens -0.0425 0.0025 3127 -16.8501 0
prev~n*time_crtl_popdens -0.0031 0.0026 3127 -1.1684 0.2428
prev~o*time_crtl_popdens*time 1e-04 0.0041 3126 0.0267 0.9787
prev~c*time_crtl_popdens*time 0.0039 0.0035 3126 1.1107 0.2668
prev~e*time_crtl_popdens*time 0.0047 0.0031 3126 1.5126 0.1305
prev~a*time_crtl_popdens*time -0.0025 0.0036 3126 -0.6789 0.4972
prev~n*time_crtl_popdens*time -0.0055 0.0024 3126 -2.2661 0.0235
# social distancing
models_socdist <- list(models_o_socdist,
models_c_socdist,
models_e_socdist,
models_a_socdist,
models_n_socdist)
sum_tab_socdist <- summary_table(models_socdist, dv_name = 'socdist', prev=T)
write.table(sum_tab_socdist, quote=F)
Value Std.Error DF t-value p-value
socdist~o*time_crtl_popdens 0.0068 0.0018 2791 3.7749 2e-04
socdist~c*time_crtl_popdens -0.006 0.0018 2791 -3.3046 0.001
socdist~e*time_crtl_popdens 0.0073 0.0018 2791 4.0335 1e-04
socdist~a*time_crtl_popdens -0.0071 0.0018 2791 -3.9228 1e-04
socdist~n*time_crtl_popdens -9e-04 0.0018 2791 -0.5053 0.6134
socdist~o*time_crtl_popdens*time -0.0012 0.003 2790 -0.4016 0.688
socdist~c*time_crtl_popdens*time 0.0015 0.0026 2790 0.548 0.5837
socdist~e*time_crtl_popdens*time 0.0026 0.0023 2790 1.1082 0.2679
socdist~a*time_crtl_popdens*time -8e-04 0.0027 2790 -0.2856 0.7752
socdist~n*time_crtl_popdens*time -0.0012 0.0018 2790 -0.6676 0.5044
socdist~o*time_crtl_popdens*time_prev -0.0016 0.0029 2789 -0.5482 0.5836
socdist~c*time_crtl_popdens*time_prev 6e-04 0.0025 2789 0.2296 0.8184
socdist~e*time_crtl_popdens*time_prev 0.0018 0.0022 2789 0.806 0.4203
socdist~a*time_crtl_popdens*time_prev -5e-04 0.0026 2789 -0.1883 0.8507
socdist~n*time_crtl_popdens*time_prev -6e-04 0.0017 2789 -0.3241 0.7459
# slope prevalence
df_uk_slope_prev <- df_uk %>% split(.$ut_area) %>%
map(~ lm(rate_day ~ time, data = .)) %>%
map(coef) %>%
map_dbl('time') %>%
as.data.frame() %>%
rownames_to_column('ut_area') %>%
rename(slope_prev = '.')
# merge with control variables
df_uk_slope_prev <- df_uk %>% select(-time, -rate_day) %>%
distinct() %>%
inner_join(df_uk_slope_prev, by = 'ut_area') %>%
drop_na()
head(df_uk_slope_prev)
NA
df_uk_slope_prev %>% ggplot(aes(slope_prev)) + geom_histogram(bins = 100)
df_uk_slope_socdist %>% ggplot(aes(slope_socdist)) + geom_histogram(bins = 100)
df_uk_slope_prev
ctrls <- cforest_unbiased(ntree=500, mtry=5)
crf_o_fit_prev <- cforest(slope_prev ~ pers_o + airport_dist + males +
popdens + manufacturing + tourism +
health + academic + medinc + medage + conservative,
df_uk_slope_prev[-1],
controls = ctrls)
crf_o_varimp_prev <- varimp(crf_o_fit_prev, nperm = 1)
crf_o_varimp_cond_prev <- varimp(crf_o_fit_prev, conditional = T, nperm = 1)
crf_o_varimp_prev %>% as.data.frame() %>% rownames_to_column('variable') %>%
ggplot(aes(x=variable, y=.)) +
geom_bar(stat = 'identity') +
theme(axis.text.x = element_text(angle = 90))
crf_o_varimp_cond_prev %>% as.data.frame() %>% rownames_to_column('variable') %>%
ggplot(aes(x=variable, y=.)) +
geom_bar(stat = 'identity') +
theme(axis.text.x = element_text(angle = 90))
ctrls <- cforest_unbiased(ntree=500, mtry=5)
crf_c_fit_prev <- cforest(slope_prev ~ pers_c + airport_dist + males +
popdens + manufacturing + tourism +
health + academic + medinc + medage + conservative,
df_uk_slope_prev[-1],
controls = ctrls)
crf_c_varimp_prev <- varimp(crf_c_fit_prev, nperm = 1)
crf_c_varimp_cond_prev <- varimp(crf_c_fit_prev, conditional = T, nperm = 1)
crf_c_varimp_prev %>% as.data.frame() %>% rownames_to_column('variable') %>%
ggplot(aes(x=variable, y=.)) +
geom_bar(stat = 'identity') +
theme(axis.text.x = element_text(angle = 90))
crf_c_varimp_cond_prev %>% as.data.frame() %>% rownames_to_column('variable') %>%
ggplot(aes(x=variable, y=.)) +
geom_bar(stat = 'identity') +
theme(axis.text.x = element_text(angle = 90))
ctrls <- cforest_unbiased(ntree=500, mtry=5)
crf_e_fit_prev <- cforest(slope_prev ~ pers_e + airport_dist + males +
popdens + manufacturing + tourism +
health + academic + medinc + medage + conservative,
df_uk_slope_prev[-1],
controls = ctrls)
crf_e_varimp_prev <- varimp(crf_e_fit_prev, nperm = 1)
crf_e_varimp_cond_prev <- varimp(crf_e_fit_prev, conditional = T, nperm = 1)
crf_e_varimp_prev %>% as.data.frame() %>% rownames_to_column('variable') %>%
ggplot(aes(x=variable, y=.)) +
geom_bar(stat = 'identity') +
theme(axis.text.x = element_text(angle = 90))
crf_e_varimp_cond_prev %>% as.data.frame() %>% rownames_to_column('variable') %>%
ggplot(aes(x=variable, y=.)) +
geom_bar(stat = 'identity') +
theme(axis.text.x = element_text(angle = 90))
ctrls <- cforest_unbiased(ntree=500, mtry=5)
crf_a_fit_prev <- cforest(slope_prev ~ pers_a + airport_dist + males +
popdens + manufacturing + tourism +
health + academic + medinc + medage + conservative,
df_uk_slope_prev[-1],
controls = ctrls)
crf_a_varimp_prev <- varimp(crf_a_fit_prev, nperm = 1)
crf_a_varimp_cond_prev <- varimp(crf_a_fit_prev, conditional = T, nperm = 1)
crf_a_varimp_prev %>% as.data.frame() %>% rownames_to_column('variable') %>%
ggplot(aes(x=variable, y=.)) +
geom_bar(stat = 'identity') +
theme(axis.text.x = element_text(angle = 90))
crf_a_varimp_cond_prev %>% as.data.frame() %>% rownames_to_column('variable') %>%
ggplot(aes(x=variable, y=.)) +
geom_bar(stat = 'identity') +
theme(axis.text.x = element_text(angle = 90))
ctrls <- cforest_unbiased(ntree=500, mtry=5)
crf_n_fit_prev <- cforest(slope_prev ~ pers_n + airport_dist + males +
popdens + manufacturing + tourism +
health + academic + medinc + medage + conservative,
df_uk_slope_prev[-1],
controls = ctrls)
crf_n_varimp_prev <- varimp(crf_n_fit_prev, nperm = 1)
crf_n_varimp_cond_prev <- varimp(crf_n_fit_prev, conditional = T, nperm = 1)
crf_n_varimp_prev %>% as.data.frame() %>% rownames_to_column('variable') %>%
ggplot(aes(x=variable, y=.)) +
geom_bar(stat = 'identity') +
theme(axis.text.x = element_text(angle = 90))
crf_n_varimp_cond_prev %>% as.data.frame() %>% rownames_to_column('variable') %>%
ggplot(aes(x=variable, y=.)) +
geom_bar(stat = 'identity') +
theme(axis.text.x = element_text(angle = 90))
Social distancing